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We report on dynamical simulations of Bose-Einstein condensation via evaporative cooling in an 
atomic trap. The results show evidence for spontaneous vortex formation and quantum dynamics 
in small traps. 
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Evaporative cooling has been successfully used to produce Bose-Einstein condensates (BEC) inside magneto-optic 
traps with neutral atoms A number of questions arise as to the quantum state that is achieved, since this involves 
Q«,^ ' both the dynamics of the cooling process and the applicability of the ergodic hypothesis. Atom-atom interactions 
^-H , have a strong influence on the cooling process and the final state in these experiments. Quantum fluctuations are 
important in determining atom laser coherence properties |^], especially since the experimental systems do not have as 
' large a particle number as traditional condensed matter experiments. However, there is no guarantee that a canonical 

■ ensemble will result from evaporative cooling, as the observations are made in a transient, non-equilibrium phase. 
] Thus, conventional canonical methods may not be applicable to these experiments. 

■ In this paper, we report the use of phase-space methods for direct quantum dynamical calculations of the cooling 
Cs| ' and formation of Bose-Einstein condensates on a three-dimensional lattice. The results are restricted as yet to small 

condensates, due to the large numbers of modes involved. The computational results are very similar to those observed 
' experimentally. In particular, we find quantum evaporative cooling, followed by a clear transition to a condensate. 
] This is strongly influenced by non-classical features of the quantum dynamics. The calculations indicate additional 
I ■ structure, which we interpret as spontaneous formation of vortices - a process of much wider interest in other fields 

of physics These appear to originate in the residual orbital angular momentum of the trapped atoms, which was 
. neglected in previous studies, and would provide a significant test of the present theory. 

' Earlier calculations of cooling dynamics have usually treated the cooling process either classically [^jj^, or have 
op used various additional assumptions about the quantum states involved. This leads to the question of how to handle 
the transition to the final quantum dominated condensate, which is often assumed to be a canonical ensemble at a 
temperature estimated from the classical theory. The final ensemble behaviour is then usually calculated from the 
mean- field Gross-Pitaevskii equations although some attempts have been made to go beyond this 0, including 
C I treatment of the kinetics of condensation based on a master equation. However, small atom traps are neither 
. in the thermodynamic limit, nor necessarily in a steady-state. A first principles theory is really needed, to provide a 
^ ' benchmark for comparisons of these previous approximations, like the quantum Monte Carlo (QMT) theories in 
O ■ equilibrium systems. 

O - In our calculations, we include 3 x 10^ relevant modes (which is a very conservative estimate), with up to 1.0 x 10* 
^ ] atoms present. The quantum state-vector therefore has over lO^'""^" components. One possible approach in principle 
is to use quantum number state calculations in the time-domain. Any direct calculation which includes all the relevant 
modes of the trapped atoms - up to the energy scales required for evaporating atoms to escape - is easily seen to be 
i an enormous computational problem. 
" " ' A more practical technique is to utilize phase-space methods that have already proved successful in laser theory. 
These techniques can handle large numbers of particles, but can also systematically treat departures from classical 
behaviour, including boson interactions. Generalized phase-space representations were used to correctly predict 
quadrature squeezed quantum soliton dynamics in optical fibers , which are described by nearly identical quantum 
equations to those used in atom-atom interaction studies. The coherent-state (positive-P) phase-space equations 
are exactly equivalent to the relevant quantum equations - provided phase-space boundary terms [ pl| vanish. They 
have the advantage that they are computationally tractable for the large Hilbert spaces typical of BEC experiments. 
Techniques of this sort can provide a first step towards extending QMT methods pO| into the time-domain. 

The model that we use includes the usual non-relativistic Hamiltonian for neutral atoms in a trap V^(x), interacting 
via a potential J7(x), together with absorbing reservoirs -R(x), in c? = 2 or d = 3 dimensions: 

[ d'^x^— V*^(x)V*(x) + F(x)v^t(x)^(x) + *t(x)i?(x) + #(x)i?t(x) + i / d'^yU (x - y)iiUx)iiUyW (yW (x) 
J l^m 2J 

Here -R(x) represents a localized absorber that removes the neutral atoms - for example, via collisions with foreign 
atoms, or at the location of the 'RF-scalpel' resonance, which is used to cause evaporative cooling pi. We expand ^ 
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using free- field modes with a momentum cut-off kmax- Provided that kmax « ^0 ^' where oq is the S-wave scattering 
length, J7(x — y) can be replaced by the renormalized pseudo-potential M(5'^(x — y), where u — 47rao?i^/TO in three 
dimensions. In two dimensions, u is defined similarly but with a factor in the denominator, which corresponds to 
the effective spatial extent of the condensate in the third direction. This factor is of the order of the lattice spacing 
in the simulation, and is chosen to be equal to a;oj the scaling length. 

The resulting quantum time-evolution for the density matrix p can be solved by expanding into a coherent-state 
basis, and then (provided phase-space boundary terms vanish) transforming to an equivalent set of equations in the 
positive-P representation. The phase-space equations in the BEC case can be expressed as two coupled complex 
partial stochastic differential equations of the form: 



dt 
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where j = 1,2 and where the stochastic fields ipj are the coherent state amplitudes of a non-diagonal coherent state 
projector, IV'i) (V-'2|/('02|V'i)- These equations can be readily simulated numerically [p^ in one, two or three transverse 
dimensions, with either attractive or repulsive potentials. The form of the potentials was chosen to be 



y(x,i) = (1 - at)F™a.Sjli[sin(^a:,/L,)]2 
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where a is typically the inverse of the total simulation time. The potential height was swept downwards linearly in 
time, thus successively removing cooler and cooler subpopulations of atoms. The absorption rate r(x) was chosen as 



r(x) = r,„,,Sli[sin(^x,/Lj)]' 



(4) 



Here Lj is the trap width in the j-th direction, such that —Lj/2 < Xj < Lj/2. The sinusoidal shape of the potential 
and absorption was chosen so that the trap would be harmonic near the center of the trap, and smoothly approach a 
maximum near the edge. Thus hot atoms are absorbed when they reach regions of large r(x), located near the trap 
edges. 

A useful feature of Eq (||) is that in the deterministic limit, this corresponds precisely to the well-known Gross- 
Pitaevskii equations, with the addition of a coefficient r(x) for the absorption of atoms by the reservoirs. Quantum 
effects come from the terms ^j, which are real Gaussian stochastic fields, with correlations: 



(ei(s,x)e2(i,y)) 
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The quantum correlations that can be calculated include n(k) = {'il)i(]i)ip2 0^)) ^ which gives the observed momentum 
distribution. 

The results of the simulations depend critically on the exact parameters chosen, just as one would expect from the 
known sensitivity of the experiments to the precise experimental conditions. In practical computations, it is necessary 
to consider rather small traps. This is because the numerical lattice spacing used to sample the stochastic fields in 
x-space must be of order Ax — I /kmax, where kmax is the largest ordinary momentum considered in the problem. 
However, the value of the corresponding kinetic energy, Ek = (7ifcmaa;)^/2m, must be large enough to allow energetic 
atoms to escape over the potential barrier of the trap, otherwise no cooling can take place. This sets an upper-bound 
on the lattice spacing, and hence on the maximum trap size - which depends on the number of lattice points that can 
be computed. 

The available lattice sizes used here were 32'' points, depending on the dimensionality d. With this limit, and 
parameter values similar to those used in the experiments, the available trap sizes that can be treated are of the order 
of micron dimensions. These are smaller than those used currently, although traps of this type are quite feasible. 
The other possibility within the constraints is to use a trap which is of larger dimensions but lower in potential 
height. For this type of trap which was simulated here, the width was Lj = lO/im, with a potential height of 
Vmax/ks — 1-9 X 10~'^K, and an initial temperature of Tq = 2.4 x 10~'^K. 

For physical reasons, a further limitation is that the initial density must be such that {n{k)) < 1 - otherwise the 
starting point would already have a Bosc-Einstcin condensation. This places a limit on the number of atoms which 
can simulated, if we assume an initially non-condensed grand canonical ensemble of (approximately) non-interacting 
atoms. There were initially around 500 atoms in the two dimensional simulations reported here, and 10, 000 in the 
three dimensional case. These corresponded to atomic densities of uq = 5.0 x 10^^/m^ and tiq — 1.0 x 10^^/m^ 
respectively. 

For the small trap parameters used in the simulations, the effect of the stochastic terms on the dynamics is very 
large. In fact the quantum fluctuations that these stochastic terms introduce are much larger than the initial thermal 
fluctuations, such that the initial features of the distribution do not persist. This means that the choice of the initial 
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state of the system is not critical, and also that in order to determine properties of the final quantum ground state of 
the system, the stochastic terms are vital. For comparison, we investigated the effect of removing the quantum noise 
terms, so that the simulations were simply of the Gross-Pitaevskii equation, with initial conditions corresponding to 
a thermal state. For our parameters, these situations did not show strong Bose condensation effects, in contrast to 
the fully quantum-mechanical simulations. This demonstrates the highly non-classical nature of the early stages of 
Bose-condensation, in which spontaneous transitions to the lowest energy states clearly play an important role. 

For the simulations shown in the figures, ao = 0.6nm and the mass, corresponding to rubidium, is m = 1.5 x lO^^^kg. 
These parameters correspond to relatively weakly interacting atoms, in order to reduce the sampling error - which 
increased rapidly with longer times and larger coupling constants. No large phase-space excursions were observed 
with these parameters. All results are plotted in normalized units, with space scaled by xq = 0.76/xm, and time scaled 
by to = 0.79ms. The time-step was typically io/2500, with all calculations being repeated at half the time-step (and 
noise sampled from the same process with twice the resolution to check numerical convergence. The boundary 

absorption term was set to Tmax — lO'^/s. 

In momentum space, the final atom density for individual trajectories in both two and three dimensions is quite 
narrow and tall, with a width corresponding to a temperature well below the critical temperature for EEC. The peak 
final momentum state population is much greater than one (and greater than the initial conditions). This is more 
pronounced in the three dimensional case than in two dimensions, showing that the evaporative cooling process is 
more efficient with the extra degree of freedom and the greater number of atoms that are present. 

As is usual in quantum mechanics, only the ensemble averages of the simulations have an operational meaning. 
Thus, while individual stochastic realizations have a definite coherent phase, these phases are different for distinct 
stochastic realizations - the ensemble average has no absolute phase information. The average evolution of a two 
dimensional condensate is shown in Fig. (1); in this case, the condensate is only weakly occupied: 

Since the condensate does not have to form in the ground-state, the Bose-condensed peaks that occur at different 
momentum values in single runs are averaged out in the overall ensemble. A more useful indication of condensation 
is given by the following measure of phase-space confinement: 

Q ^ /d^fc(V'i(k)V;i(k)v>*(k)v>,*(k)) 

This higher order correlation function is the quantum analogue of the participation ratio defined by Hall . Figure 
(2) shows the evolution of Q calculated from 15 runs of the three dimensional simulation. The sharp rise near t = 100 
is a strong indication of condensation occurring at this point. 

For the finite-size condensates in atom traps, just as the final ground state is not expected to be precisely the 
zero momentum eigenstate, so too such condensates are not constrained to fall into the J = angular momentum 
eigenstate. Both the initial and the escaping atoms have an arbitrary angular momentum. We can estimate that 
the variance in angular momentum will scale approximately as (J^) cx iV, from central limit theorem arguments. 
Thus, we can expect that each trapped condensate should have angular momentum, unless constrained by the trap 
geometry. The angular momentum can be carried either by quasi-particles or vortices, although a volume-filling j-th 
order vortex has J — Nj^ and therefore cannot form spontaneously in the thermodynamic limit of large N. For small 
condensates, a, j = ±1 vortex may be quite likely. Several authors [ p^ have considered how such vortex states may 
form through stirring or rotating a condensate, and the stability of vortices has been explored [|l5|. Here we consider 
the possibility of vortices forming spontaneously in the condensate through the process of evaporative cooling, without 
external intervention. 

The presence of vortex states can be detected quantitatively by transforming the spatial lattice into a lattice which 
uses the angular momentum eigenstates as a basis set. The two-dimensional results, which are presented here, are 
obtained by integrating the spatial profile over orthogonal modes with corresponding field operators '^jn- The angular 
momentum distribution is then given by a summation over the radial modes: 

n 

The angular momentum distribution for individual trajectories shows large occupation in particular angular modes, 
different for each run. This indicates that vortices with different momenta appear each time. For example, in one 
run, a vortex with j = — 1 appears at about one quarter of the way through the simulation, and persists until the end. 
The maximum occupation of the vortex is around n(j) = 20, owing to relatively small initial atom numbers in this 
2D trap simulation. Shown in Fig. (3) is the ensemble average of the angular momentum distribution, which reveals 
quite a broad range of final angular momentum. This is consistent with the existence of vortices. 

In summary, we have demonstrated a three-dimensional real-time quantum dynamical simulation of Bose- 
condensation with mesoscopic numbers of interacting atoms on a large lattice. Sampling errors and lattice size 
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restrictions impose strong limitations on these initial simulations. The results, as well as showing evidence for highly 
non-classical behaviour in a first principles simulation of BEC formation, indicate the possibility of spontaneous vortex 
formation in small evaporatively cooled condensates. 
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FIG. 1. Simulation of a two-dimensional Bose condensate, showing the ensemble average (55 paths) atom density {n{k)) 
along one dimension in fourier space versus time. 

FIG. 2. Simulation of a three-dimensional Bose condensate, showing the ensemble average evolution (15 paths) of the 
confinement parameter Q. 

FIG. 3. Ensemble average of the angular momentum distribution {n{j)), during the condensation of a two-dimensional 
Bose-condensate (40 paths). 
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